#
# Copyright (C) 2023, Inria
# GRAPHDECO research group, https://team.inria.fr/graphdeco
# All rights reserved.
#
# This software is free for non-commercial, research and evaluation use 
# under the terms of the LICENSE.md file.
#
# For inquiries contact  george.drettakis@inria.fr
#

import torch
import numpy as np
from utils.general_utils import inverse_sigmoid, get_expon_lr_func, build_rotation
from torch import nn
import os, math
from utils.system_utils import mkdir_p
from plyfile import PlyData, PlyElement
from utils.sh_utils import RGB2SH
from simple_knn._C import distCUDA2
from utils.graphics_utils import BasicPointCloud, hemisphere
from utils.general_utils import strip_symmetric, build_scaling_rotation
from scene.embedding import Embedding, MLP, PosEmbedding

class GaussianModel:

    def setup_functions(self):
        def build_covariance_from_scaling_rotation(scaling, scaling_modifier, rotation):
            L = build_scaling_rotation(scaling_modifier * scaling, rotation)
            actual_covariance = L @ L.transpose(1, 2)
            symm = strip_symmetric(actual_covariance)
            return symm
        
        # 尺度矩阵：S, 旋转矩阵：R
        self.scaling_activation = torch.exp
        self.scaling_inverse_activation = torch.log
        # 协方差矩阵：RS(S^T)(R^T)
        self.covariance_activation = build_covariance_from_scaling_rotation

        self.opacity_activation = torch.sigmoid
        self.inverse_opacity_activation = inverse_sigmoid

        self.rotation_activation = torch.nn.functional.normalize

    def __init__(self, args):
        self.active_sh_degree = 0               # 球谐阶数
        # 存储不同信息的张量（tensor）
        self._xyz = torch.empty(0)              # 空间位置
        self._features = torch.empty(0)         # 特征属性
        self._scaling = torch.empty(0)          # 椭球的形状尺度
        self._rotation = torch.empty(0)         # 椭球的旋转方向
        self._opacity = torch.empty(0)          # 椭球的透明度
        
        self.max_radii2D = torch.empty(0)       # 投影到相机平面的椭圆的最大包围半径
        
        self.xyz_gradient_accum = torch.empty(0)    # 空间位置的累计梯度
        self.denom = torch.empty(0)                 # 记录每个Gaussian被用到的次数，为了计算平均梯度
                
        self.optimizer = None                   # 初始化的优化器为 None
        self.percent_dense = 0                  # 初始化密度的的百分比为 0
        self.spatial_lr_scale = 0               # 初始化的空间学习率缩放为 0
        self.setup_functions()                  # 初始化不同信息的张量的激活函数
        
        self.vis_emb_dim = args.vis_emb_dim
        self.dir_emb_dim = args.dir_emb_dim
        self.sun_emb_dim = args.sun_emb_dim
        self.env_emb_dim = args.env_emb_dim
        
        self.vis_feat_dim = args.vis_feat_dim
        self.ref_feat_dim = args.ref_feat_dim
        self.sun_feat_dim = args.sun_feat_dim
        self.env_feat_dim = args.env_feat_dim
        
        self.feat_dim = self.ref_feat_dim + self.vis_feat_dim + self.sun_feat_dim + self.env_feat_dim
        self.hidden_dim = args.hidden_dim
        
        ref_dim_to_mlp = self.ref_feat_dim
        vis_dim_to_mlp = self.vis_feat_dim + self.vis_emb_dim
        sun_dim_to_mlp = self.sun_feat_dim + self.sun_emb_dim
        env_dim_to_mlp = self.env_feat_dim + self.env_emb_dim
        
        self._ref_mlp = MLP(ref_dim_to_mlp, self.hidden_dim, 3, 2).cuda()
        self._vis_mlp = MLP(vis_dim_to_mlp, self.hidden_dim, 1, 2).cuda()
        self._sun_mlp = MLP(sun_dim_to_mlp, self.hidden_dim, 3, 2).cuda()
        self._env_mlp = MLP(env_dim_to_mlp, self.hidden_dim, 3, 2).cuda()
        
                            
    def feature_encoder(self, camera, visible_mask=None, vis_idx=None, sun_idx=None, env_idx=None):
        if visible_mask is None:
            visible_mask = torch.full(self.get_xyz.shape[0], True).cuda()
        if vis_idx is None:
            vis_idx = camera.uid
        if sun_idx is None:
            sun_idx = camera.uid
        if env_idx is None:
            env_idx = camera.uid
            
        feats = self.get_features[visible_mask]
        
        ref_feat, vis_feat, sun_feat, env_feat = \
            torch.split(feats, (self.ref_feat_dim, self.vis_feat_dim, self.sun_feat_dim, self.env_feat_dim), dim=-1)
        # ref_feat, sun_feat, env_feat = torch.split(feats, (self.ref_feat_dim, self.sun_feat_dim, self.env_feat_dim), dim=-1)
        
        mlps = self.get_mlps
        embeddings = self.get_embeddings
        tem_idx = torch.ones_like(feats[:, 0], dtype=torch.long)
        
        # reflectance
        reflectance = mlps["ref_mlp"](ref_feat)
        
        # sun visibility
        if isinstance(vis_idx, dict):
            idx0, idx1, weight = tem_idx*vis_idx['idx0'], tem_idx*vis_idx['idx1'], vis_idx['weight']
            vis_emb = embeddings["vis_embedding"](idx0) * (1-weight) + embeddings["vis_embedding"](idx1) * weight
        else:
            vis_idx = tem_idx * vis_idx
            vis_emb = embeddings["vis_embedding"](vis_idx)   # (N, 3)
        vis_feat = torch.cat([vis_feat, vis_emb], -1)
        visibility = mlps["vis_mlp"](vis_feat)
        
        # sun shade
        if isinstance(sun_idx, dict):
            idx0, idx1, weight = tem_idx*sun_idx['idx0'], tem_idx*sun_idx['idx1'], sun_idx['weight']
            sun_emb = embeddings["sun_embedding"](idx0) * (1-weight) + embeddings["sun_embedding"](idx1) * weight
        else:
            sun_idx = tem_idx * sun_idx
            sun_emb = embeddings["sun_embedding"](sun_idx)   # (N, sun_feat_dim)
        sun_feat = torch.cat([sun_feat, sun_emb], -1)
        shade_sun = mlps["sun_mlp"](sun_feat)
        
        # env shade
        if isinstance(env_idx, dict):
            idx0, idx1, weight = tem_idx*env_idx['idx0'], tem_idx*env_idx['idx1'], env_idx['weight']
            env_emb = embeddings["env_embedding"](idx0) * (1-weight) + embeddings["env_embedding"](idx1) * weight
        else:
            env_idx = tem_idx * env_idx
            env_emb = embeddings["env_embedding"](env_idx)   # (N, sun_feat_dim)
        env_feat = torch.cat([env_feat, env_emb], -1)
        shade_env = mlps["env_mlp"](env_feat)
        
        del ref_feat, vis_feat, sun_feat, env_feat, sun_emb, env_emb
        return torch.cat([reflectance, visibility, shade_sun, shade_env], dim=-1)
                    
    def init_embeddings(self, num_cameras):
        self._vis_embedding = torch.nn.Embedding(num_cameras, self.vis_emb_dim)
        self._vis_embedding.weight = nn.Parameter(torch.randn(num_cameras, self.vis_emb_dim) / math.sqrt(self.vis_emb_dim/2))
        self._vis_embedding = self._vis_embedding.cuda()
        
        # self._dir_embedding = torch.nn.Embedding(num_cameras, self.dir_emb_dim)
        # self._dir_embedding.weight = nn.Parameter(hemisphere(torch.randn(num_cameras, self.dir_emb_dim)))
        # self._dir_embedding = self._dir_embedding.cuda()
        
        self._sun_embedding = torch.nn.Embedding(num_cameras, self.sun_emb_dim)
        self._sun_embedding.weight = nn.Parameter(torch.randn(num_cameras, self.sun_emb_dim) / math.sqrt(self.sun_emb_dim/2))
        self._sun_embedding = self._sun_embedding.cuda()
        
        self._env_embedding = torch.nn.Embedding(num_cameras, self.env_emb_dim)
        self._env_embedding.weight = nn.Parameter(torch.randn(num_cameras, self.env_emb_dim) / math.sqrt(self.env_emb_dim/2))
        self._env_embedding = self._env_embedding.cuda()        
    
    # 获取 3D Gaussian 的属性信息和优化器信息
    def capture(self):
        return (
            self.active_sh_degree,
            self._xyz,
            self._features,
            self._scaling,
            self._rotation,
            self._opacity,
            self.max_radii2D,
            self.xyz_gradient_accum,
            self.denom,
            self.optimizer.state_dict(),
            self.spatial_lr_scale,
        )
    
    # 保存训练信息和优化后的 3DGS 的属性信息
    def restore(self, model_args, training_args):
        (self.active_sh_degree, 
        self._xyz, 
        self._features, 
        self._scaling, 
        self._rotation, 
        self._opacity,
        self.max_radii2D, 
        xyz_gradient_accum, 
        denom,
        opt_dict, 
        self.spatial_lr_scale) = model_args
        self.training_setup(training_args)
        self.xyz_gradient_accum = xyz_gradient_accum
        self.denom = denom
        self.optimizer.load_state_dict(opt_dict)

    # 获取激活后的椭球的形状尺度
    @property
    def get_scaling(self):
        return self.scaling_activation(self._scaling)
    
    # 获取激活后的椭球的旋转方向
    @property
    def get_rotation(self):
        return self.rotation_activation(self._rotation)
    
    # 获取椭球的空间位置
    @property
    def get_xyz(self):
        return self._xyz
    
    # 获取椭球的各阶球谐系数
    @property
    def get_features(self):
        return self._features
    
    # 获取激活后的椭球的透明度
    @property
    def get_opacity(self):
        return self.opacity_activation(self._opacity)
    
    @property
    def get_mlps(self):
        return {'ref_mlp': self._ref_mlp,
                'vis_mlp': self._vis_mlp,
                'sun_mlp': self._sun_mlp,
                'env_mlp': self._env_mlp}
        
    @property
    def get_embeddings(self):
        return {'vis_embedding': self._vis_embedding,
                # 'dir_embedding': self._dir_embedding,
                'sun_embedding': self._sun_embedding,
                'env_embedding': self._env_embedding}
            
    # 获取激活后的椭球的协方差矩阵
    def get_covariance(self, scaling_modifier=1):
        return self.covariance_activation(self.get_scaling, scaling_modifier, self._rotation)

    # 经过一定轮次的迭代之后，球谐阶数 +1 （最大球谐阶数 = self.max_sh_degree）
    def oneupSHdegree(self):
        if self.active_sh_degree < self.max_sh_degree:
            self.active_sh_degree += 1

    # 从给定或者随机初始化的点云数据中初始化 3D Gaussians
    def create_from_pcd(self, pcd : BasicPointCloud, spatial_lr_scale : float):
        self.spatial_lr_scale = spatial_lr_scale
        
        # 从点云数据中构建空间点的张量：(N, 3)
        fused_point_cloud = torch.tensor(np.asarray(pcd.points)).float().cuda()
        # 初始化可学习特征
        features = (torch.randn(fused_point_cloud.shape[0], self.feat_dim) / math.sqrt(self.feat_dim/2)).cuda()

        print("Number of points at initialisation : ", fused_point_cloud.shape[0])

        # distCUDA2 计算点云中的每个点到与其最近的 K 个点的平均距离的平方：(N, )
        dist2 = torch.clamp_min(distCUDA2(torch.from_numpy(np.asarray(pcd.points)).float().cuda()), 0.0000001)
        # 每个点在 x, y, z 方向上的尺度：(N, 3)
        scales = torch.log(torch.sqrt(dist2))[...,None].repeat(1, 3)        # (N, 3)
        # 每个点的旋转参数，四元数组：(N, 4)
        rots = torch.zeros((fused_point_cloud.shape[0], 4), device="cuda")  # (N, 4)
        rots[:, 0] = 1

        # 每个 Gaussian 的不透明度，初始化为 0.1：(N, 1)
        opacities = inverse_sigmoid(0.1 * torch.ones((fused_point_cloud.shape[0], 1), dtype=torch.float, device="cuda"))

        self._xyz = nn.Parameter(fused_point_cloud.requires_grad_(True))    # (N, 3)
        self._features = nn.Parameter(features.requires_grad_(True))        # (N, 48)
        self._scaling = nn.Parameter(scales.requires_grad_(True))           # (N, 3)
        self._rotation = nn.Parameter(rots.requires_grad_(True))            # (N, 4)
        self._opacity = nn.Parameter(opacities.requires_grad_(True))        # (N, 1)
        self.max_radii2D = torch.zeros((self.get_xyz.shape[0]), device="cuda")   # (N, )

    def training_setup(self, training_args):
        self.percent_dense = training_args.percent_dense # 0.01
        # 存储每个 3D Gaussian 的均值 xyz 的梯度，用于判断是否对 3D Gaussian 进行稠密化
        self.xyz_gradient_accum = torch.zeros((self.get_xyz.shape[0], 1), device="cuda")    # (N, 1)
        self.denom = torch.zeros((self.get_xyz.shape[0], 1), device="cuda")                 # (N, 1)

        l = [
            {'params': [self._xyz], 'lr': training_args.position_lr_init * self.spatial_lr_scale, "name": "xyz"},
            {'params': [self._features], 'lr': training_args.feature_lr, "name": "features"},
            {'params': [self._opacity], 'lr': training_args.opacity_lr, "name": "opacity"},
            {'params': [self._scaling], 'lr': training_args.scaling_lr, "name": "scaling"},
            {'params': [self._rotation], 'lr': training_args.rotation_lr, "name": "rotation"},
            
            {'params': self._ref_mlp.parameters(), 'lr': training_args.mlp_lr_init, "name": "ref_mlp"},
            {'params': self._vis_mlp.parameters(), 'lr': training_args.mlp_lr_init, "name": "vis_mlp"},
            {'params': self._env_mlp.parameters(), 'lr': training_args.mlp_lr_init, "name": "env_mlp"},
            {'params': self._sun_mlp.parameters(), 'lr': training_args.mlp_lr_init, "name": "sun_mlp"},
            
            {'params': self._vis_embedding.parameters(), 'lr': training_args.embedding_lr_init, "name": "vis_embedding"},
            # {'params': self._dir_embedding.parameters(), 'lr': training_args.embedding_lr_init, "name": "dir_embedding"},
            {'params': self._sun_embedding.parameters(), 'lr': training_args.embedding_lr_init, "name": "sun_embedding"},
            {'params': self._env_embedding.parameters(), 'lr': training_args.embedding_lr_init, "name": "env_embedding"}
        ]

        self.optimizer = torch.optim.Adam(l, lr=0.0, eps=1e-15)
        self.xyz_scheduler_args = get_expon_lr_func(lr_init=training_args.position_lr_init*self.spatial_lr_scale,
                                                    lr_final=training_args.position_lr_final*self.spatial_lr_scale,
                                                    lr_delay_mult=training_args.position_lr_delay_mult,
                                                    max_steps=training_args.position_lr_max_steps)
        
        self.mlp_scheduler_args = get_expon_lr_func(lr_init=training_args.mlp_lr_init,
                                                    lr_final=training_args.mlp_lr_final,
                                                    lr_delay_mult=training_args.mlp_lr_delay_mult,
                                                    max_steps=training_args.mlp_lr_max_steps)
        
        self.emb_scheduler_args = get_expon_lr_func(lr_init=training_args.embedding_lr_init,
                                                    lr_final=training_args.embedding_lr_final,
                                                    lr_delay_mult=training_args.embedding_lr_delay_mult,
                                                    max_steps=training_args.embedding_lr_max_steps)
        
    def update_learning_rate(self, iteration):
        ''' Learning rate scheduling per step '''
        for param_group in self.optimizer.param_groups:
            if param_group["name"] == "xyz":
                lr = self.xyz_scheduler_args(iteration)
                param_group['lr'] = lr
            elif param_group["name"].endswith("mlp"):
                lr = self.mlp_scheduler_args(iteration)
                param_group['lr'] = lr
            elif param_group["name"].endswith("embedding"):
                lr = self.emb_scheduler_args(iteration)
                param_group['lr'] = lr

    def construct_list_of_attributes(self):
        l = ['x', 'y', 'z', 'nx', 'ny', 'nz']
        # All channels except the 3 DC
        for i in range(self._features.shape[1]):
            l.append('feat_{}'.format(i))
        l.append('opacity')
        for i in range(self._scaling.shape[1]):
            l.append('scale_{}'.format(i))
        for i in range(self._rotation.shape[1]):
            l.append('rot_{}'.format(i))
        return l

    def save_ply(self, path):
        mkdir_p(os.path.dirname(path))

        xyz = self._xyz.detach().cpu().numpy()  # (N, 3)
        normals = np.zeros_like(xyz)            # (N, 3)
        feat = self._features.detach().cpu().numpy()        # (N, 48)
        opacities = self._opacity.detach().cpu().numpy()    # (N, 1)
        scale = self._scaling.detach().cpu().numpy()        # (N, 3)
        rotation = self._rotation.detach().cpu().numpy()    # (N, 4)

        dtype_full = [(attribute, 'f4') for attribute in self.construct_list_of_attributes()]

        elements = np.empty(xyz.shape[0], dtype=dtype_full)
        attributes = np.concatenate((xyz, normals, feat, opacities, scale, rotation), axis=1)
        elements[:] = list(map(tuple, attributes))
        el = PlyElement.describe(elements, 'vertex')
        PlyData([el]).write(path)

    # 重置不透明度
    def reset_opacity(self):
        opacities_new = inverse_sigmoid(torch.min(self.get_opacity, torch.ones_like(self.get_opacity)*0.01))
        optimizable_tensors = self.replace_tensor_to_optimizer(opacities_new, "opacity")
        self._opacity = optimizable_tensors["opacity"]

    def load_ply(self, path):
        plydata = PlyData.read(path)

        xyz = np.stack((np.asarray(plydata.elements[0]["x"]),
                        np.asarray(plydata.elements[0]["y"]),
                        np.asarray(plydata.elements[0]["z"])),  axis=1)
        opacities = np.asarray(plydata.elements[0]["opacity"])[..., np.newaxis]
            
        f_names = [p.name for p in plydata.elements[0].properties if p.name.startswith("feat_")]
        f_names = sorted(f_names, key = lambda x: int(x.split('_')[-1]))
        features = np.zeros((xyz.shape[0], self.feat_dim))
        for idx, attr_name in enumerate(f_names):
            features[:, idx] = np.asarray(plydata.elements[0][attr_name])

        scale_names = [p.name for p in plydata.elements[0].properties if p.name.startswith("scale_")]
        scale_names = sorted(scale_names, key = lambda x: int(x.split('_')[-1]))
        scales = np.zeros((xyz.shape[0], len(scale_names)))
        for idx, attr_name in enumerate(scale_names):
            scales[:, idx] = np.asarray(plydata.elements[0][attr_name])

        rot_names = [p.name for p in plydata.elements[0].properties if p.name.startswith("rot")]
        rot_names = sorted(rot_names, key = lambda x: int(x.split('_')[-1]))
        rots = np.zeros((xyz.shape[0], len(rot_names)))
        for idx, attr_name in enumerate(rot_names):
            rots[:, idx] = np.asarray(plydata.elements[0][attr_name])

        self._xyz = nn.Parameter(torch.tensor(xyz, dtype=torch.float, device="cuda").requires_grad_(True))
        self._features = nn.Parameter(torch.tensor(features, dtype=torch.float, device="cuda").requires_grad_(True))
        self._opacity = nn.Parameter(torch.tensor(opacities, dtype=torch.float, device="cuda").requires_grad_(True))
        self._scaling = nn.Parameter(torch.tensor(scales, dtype=torch.float, device="cuda").requires_grad_(True))
        self._rotation = nn.Parameter(torch.tensor(rots, dtype=torch.float, device="cuda").requires_grad_(True))

    def replace_tensor_to_optimizer(self, tensor, name):
        optimizable_tensors = {}
        for group in self.optimizer.param_groups:
            if group["name"].endswith('mlp') or group["name"].endswith('embedding'):
                continue
            if group["name"] == name:
                stored_state = self.optimizer.state.get(group['params'][0], None)
                stored_state["exp_avg"] = torch.zeros_like(tensor)
                stored_state["exp_avg_sq"] = torch.zeros_like(tensor)

                del self.optimizer.state[group['params'][0]]
                group["params"][0] = nn.Parameter(tensor.requires_grad_(True))
                self.optimizer.state[group['params'][0]] = stored_state

                optimizable_tensors[group["name"]] = group["params"][0]
        return optimizable_tensors

    def _prune_optimizer(self, mask):
        optimizable_tensors = {}
        for group in self.optimizer.param_groups:
            if group["name"].endswith('mlp') or group["name"].endswith('embedding'):
                continue
            stored_state = self.optimizer.state.get(group['params'][0], None)
            if stored_state is not None:
                stored_state["exp_avg"] = stored_state["exp_avg"][mask]
                stored_state["exp_avg_sq"] = stored_state["exp_avg_sq"][mask]

                del self.optimizer.state[group['params'][0]]
                group["params"][0] = nn.Parameter((group["params"][0][mask].requires_grad_(True)))
                self.optimizer.state[group['params'][0]] = stored_state

                optimizable_tensors[group["name"]] = group["params"][0]
            else:
                group["params"][0] = nn.Parameter(group["params"][0][mask].requires_grad_(True))
                optimizable_tensors[group["name"]] = group["params"][0]
        return optimizable_tensors

    # 删除不符合要求的3D gaussian 在self.optimizer 中对应的参数(均值、球谐系数、不透明度、尺度、旋转参数)
    def prune_points(self, mask):
        valid_points_mask = ~mask
        optimizable_tensors = self._prune_optimizer(valid_points_mask)

        self._xyz = optimizable_tensors["xyz"]
        self._features = optimizable_tensors["features"]
        self._opacity = optimizable_tensors["opacity"]
        self._scaling = optimizable_tensors["scaling"]
        self._rotation = optimizable_tensors["rotation"]

        self.xyz_gradient_accum = self.xyz_gradient_accum[valid_points_mask]
        self.denom = self.denom[valid_points_mask]
        self.max_radii2D = self.max_radii2D[valid_points_mask]

    def cat_tensors_to_optimizer(self, tensors_dict):
        optimizable_tensors = {}
        for group in self.optimizer.param_groups:
            if group["name"].endswith('mlp') or group["name"].endswith('embedding'):
                continue
            assert len(group["params"]) == 1
            extension_tensor = tensors_dict[group["name"]]
            stored_state = self.optimizer.state.get(group['params'][0], None)
            if stored_state is not None:

                stored_state["exp_avg"] = torch.cat((stored_state["exp_avg"], torch.zeros_like(extension_tensor)), dim=0)
                stored_state["exp_avg_sq"] = torch.cat((stored_state["exp_avg_sq"], torch.zeros_like(extension_tensor)), dim=0)

                del self.optimizer.state[group['params'][0]]
                group["params"][0] = nn.Parameter(torch.cat((group["params"][0], extension_tensor), dim=0).requires_grad_(True))
                self.optimizer.state[group['params'][0]] = stored_state

                optimizable_tensors[group["name"]] = group["params"][0]
            else:
                group["params"][0] = nn.Parameter(torch.cat((group["params"][0], extension_tensor), dim=0).requires_grad_(True))
                optimizable_tensors[group["name"]] = group["params"][0]

        return optimizable_tensors

    # 将需要进行 clone 的 3D Gaussian 属性进行 clone，并与原 3D Gaussian 进行合并
    def densification_postfix(self, new_xyz, new_features, new_opacities, new_scaling, new_rotation):
        # 需要进行 clone 的 3D Gaussian 属性
        d = {"xyz": new_xyz,
            "features": new_features,
            "opacity": new_opacities,
            "scaling" : new_scaling,
            "rotation" : new_rotation}

        # 将属性进行 clone 并合并进行优化
        optimizable_tensors = self.cat_tensors_to_optimizer(d)
        self._xyz = optimizable_tensors["xyz"]              # (N+P, 3)
        self._features = optimizable_tensors["features"]        # (N+P, 48)
        self._opacity = optimizable_tensors["opacity"]      # (N+P, 1)
        self._scaling = optimizable_tensors["scaling"]      # (N+P, 3)
        self._rotation = optimizable_tensors["rotation"]    # (N+P, 4)

        # clone 之后重新初始化对应的均值梯度。。。
        self.xyz_gradient_accum = torch.zeros((self.get_xyz.shape[0], 1), device="cuda")
        self.denom = torch.zeros((self.get_xyz.shape[0], 1), device="cuda")
        self.max_radii2D = torch.zeros((self.get_xyz.shape[0]), device="cuda")
 
    def densify_and_split(self, grads, grad_threshold, scene_extent, N=2):
        n_init_points = self.get_xyz.shape[0]
        padded_grad = torch.zeros((n_init_points), device="cuda")
        padded_grad[:grads.shape[0]] = grads.squeeze()
        # 筛选条件：均值的梯度过大 and 尺度过大
        # 1、提取满足梯度条件的点：均值的梯度大于设定的梯度阈值
        selected_pts_mask = torch.where(padded_grad >= grad_threshold, True, False)
        # 2、提取满足尺度条件的点：尺度大于设定的尺度阈值
        selected_pts_mask = torch.logical_and(selected_pts_mask,
                                              torch.max(self.get_scaling, dim=1).values > self.percent_dense*scene_extent)
        # 把筛选出来的尺度 tensor 复制 N 次
        stds = self.get_scaling[selected_pts_mask].repeat(N,1)  # (P, 3)
        means = torch.zeros((stds.size(0), 3),device="cuda")    # (P, 3)
        # 构建均值为 0，方差为满足筛选条件的点的尺度值的新的点
        samples = torch.normal(mean=means, std=stds)            # (P, 3)
        # 将 N x 4 的旋转四元组转换成 N x 3 x 3 的旋转矩阵, N 为点云中点的个数或者3D gaussian的个数
        rots = build_rotation(self._rotation[selected_pts_mask]).repeat(N,1,1)
        # 在以原来 3D Gaussian 的均值 xyz 为中心, stds 为形状, rots 为方向的椭球内随机采样新的 3D Gaussian
        new_xyz = torch.bmm(rots, samples.unsqueeze(-1)).squeeze(-1) + self.get_xyz[selected_pts_mask].repeat(N, 1)
        # 尺度属性进行 split 之后，将每个尺度适当缩小
        new_scaling = self.scaling_inverse_activation(self.get_scaling[selected_pts_mask].repeat(N,1) / (0.8*N))
        new_rotation = self._rotation[selected_pts_mask].repeat(N,1)
        new_features = self._features[selected_pts_mask].repeat(N,1)
        new_opacity = self._opacity[selected_pts_mask].repeat(N,1)

        self.densification_postfix(new_xyz, new_features, new_opacity, new_scaling, new_rotation)

        # 将原来的那些均值的梯度超过一定阈值且尺度大于一定阈值的3D gaussian进行删除 
        # (因为已经将它们分割成了两个新的3D gaussian，原先的不再需要了)
        prune_filter = torch.cat((selected_pts_mask, torch.zeros(N * selected_pts_mask.sum(), device="cuda", dtype=bool)))
        self.prune_points(prune_filter)

    def densify_and_clone(self, grads, grad_threshold, scene_extent):
        # 筛选条件：均值的梯度过大 and 尺度过小
        # 1、提取满足梯度条件的点：均值的梯度大于设定的梯度阈值
        selected_pts_mask = torch.where(torch.norm(grads, dim=-1) >= grad_threshold, True, False)
        # 2、提取满足尺度条件的点：尺度小于设定的尺度阈值
        selected_pts_mask = torch.logical_and(selected_pts_mask,
                                              torch.max(self.get_scaling, dim=1).values <= self.percent_dense*scene_extent)
        # 筛选出满足条件的 3D Gaussian 的属性特征
        new_xyz = self._xyz[selected_pts_mask]                      # (P, 3)
        new_features = self._features[selected_pts_mask]            # (P, 48)
        new_opacities = self._opacity[selected_pts_mask]            # (P, 1)
        new_scaling = self._scaling[selected_pts_mask]              # (P, 3)
        new_rotation = self._rotation[selected_pts_mask]            # (P, 4)

        self.densification_postfix(new_xyz, new_features, new_opacities, new_scaling, new_rotation)

    # 根据梯度对 3D Gaussian 进行增加或删除
    def densify_and_prune(self, max_grad, min_opacity, extent, max_screen_size):
        xyz_grads = self.xyz_gradient_accum / self.denom    # 3D Gaussian 的均值的累积梯度
        xyz_grads[xyz_grads.isnan()] = 0.0  # 修正为 NaN 的梯度值为 0
                
        # 如果某些 3D Gaussian 的均值的梯度过大且尺度小于一定阈值，说明是欠重建，则将其进行克隆
        self.densify_and_clone(xyz_grads, max_grad, extent)
        # 如果某些 3D Gaussian 的均值的梯度过小且尺度大于一定阈值，说明是过重建，则将其进行分割
        self.densify_and_split(xyz_grads, max_grad, extent)

        # 删除不透明度小于设定阈值 and 椭球体半径过大的 3D Gaussian
        # 1、选择不透明度小于一定阈值的 3D Gaussian
        prune_mask = (self.get_opacity < min_opacity).squeeze()
        # 2、选择椭球体半径过大的 3D Gaussian
        if max_screen_size:
            big_points_vs = self.max_radii2D > max_screen_size
            big_points_ws = self.get_scaling.max(dim=1).values > 0.1 * extent
            prune_mask = torch.logical_or(torch.logical_or(prune_mask, big_points_vs), big_points_ws)
        self.prune_points(prune_mask)

        torch.cuda.empty_cache()

    def add_densification_stats(self, viewspace_point_tensor, update_filter):
        # 记录视锥内的 xyz 的梯度
        self.xyz_gradient_accum[update_filter] += torch.norm(viewspace_point_tensor.grad[update_filter,:2], dim=-1, keepdim=True)
        # 统计每个 3D Gaussian 的 xyz 的梯度被更新的次数，后续进行梯度球平均
        self.denom[update_filter] += 1
        
    def save_mlp_checkpoints(self, path):
        mkdir_p(os.path.dirname(path))
        mlps = self.get_mlps
        embeddings = self.get_embeddings   
        ckpt = {k: v.state_dict() for k, v in mlps.items()}
        ckpt.update({k: v.state_dict() for k, v in embeddings.items()})
        torch.save(ckpt, os.path.join(path, 'checkpoints.pth'))

    def load_mlp_checkpoints(self, path):
        ckpt = torch.load(os.path.join(path, 'checkpoints.pth'))
        
        mlps = self.get_mlps
        embeddings = self.get_embeddings
        for key in mlps.keys():
            mlps[key].load_state_dict(ckpt[key])
        for key in embeddings.keys():
            embeddings[key].load_state_dict(ckpt[key])
        
